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Chapter  1 


Introduction 

1.1  Rotor  Analysis  in  Ground  Effect 

The  wake  shed  by  a  helicopter  rotor  is  one  of  the  most  important  factors  in  determin¬ 
ing  rotor  performance.  Unfortunately,  the  rotor  wake  is  an  extremely  complicated 
aerodynamic  feature.  The  saying  that  “a  single  helicopter  blade  faces  more  aero¬ 
dynamic  problems  in  one  revolution  than  a  fixed  wing  aircraft  meets  in  its  entire 
lifetime”  [1]  is  particularly  true  when  the  blade  in  question  is  in  proximity  to  the 
ground,  that  is,  when  the  rotor  is  “in  ground  effect”.  In  ground  effect,  the  rotor  wake 
becomes  more  complicated.  A  low  advance  ratio  (ratio  of  free-stream  or  cross-wind  to 
rotor  tip  speed)  makes  the  rotor  wake  problem  more  complicated  still  while  retaining 
the  first-order  impact  of  wake  geometry  on  performance.  Under  these  conditions, 
the  free-stream  and/or  cross-stream  interacts  with  the  already  complex  rotor  wake  to 
produce  a  powerful  ground  vortex.  This  ground  vortex  (visualized  in  Figure  1.1  (6) 
using  helium  bubbles)  wraps  around  the  helicopter  and  causes  an  array  of  problems 
including  loss  of  yaw  control.  -In  other,  cases,  the  ground  vortex  can  interact  negative- 
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1.1.  ROTOR  ANALYSIS  IN  GROUND  EFFECT 


Figure  1.1:  Ground  Vortex  Impinging  on  Empennage 


ly  with  the  main  rotor  or  with  other  aircraft  in  formation.  Also,  the  helicopter  wake 
can  interact  with  ground  operations  or  ship  deck  operations. 

In  general,  the  computation  of  rotors  in  near-ground  operation  requires  the  ability 
to  predict  both  the  wake  structures  shed  by  the  blades  and  the  larger  vortex  structures 
along  the  ground.  It  is  the  purpose  of  this  paper  to  investigate  the  effectiveness  of 
the  Vorticity  Embedding  method  coupled  with  a  conventional  Euler  solver  to  predict 
the  wake  and  ground  vortex  structures  mentioned  above.  An  existing  but  previously 
incomplete  and  undocumented  code  (called  Helix-III)  will  be  used  for  this  purpose. 


2 


1.2.  REVIEW  OF  EXISTING  METHODS 


Figure  1.2:  Smoke  Visualization  of  Rotor  Wake  (taken  from  [11]) 


1.2  Review  of  Existing  Methods 


The  flow  field  created  by  a  rotor  is  characterized  by  thin  regions  of  vorticity  in  an 
otherwise  potential  flow.  The  smoke  visualization  of  a  rotor  wake  in  Figure  1.2  which 
is  diagrammed  in  Figure  1.3  clearly  shows  the  velocity  discontinuity  caused  by  the 
thin  region  of  vorticity  which  comprises  the  wake.  It  also  clearly  shows  the  powerful 
tip  vortices  caused  by  the  roll-up  of  the  rotor  wake.  Many  methods  have  been  devised 
to  treat  the  problems  which  are  associated  with  calculating  flows  containing  regions 
of  compact  vorticity. 


1-2.1  Vortex  Methods 

Vortex  methods  treat  the  wake  as  a  sheet  of  Lagrangian  vortex  filaments  using  the 
Biot-Savart  law  to  calculate  the  velocity  induced  by  the  filaments.  Determining  ve- 
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1.2.  REVIEW  OF  EXISTING  METHODS 


Figure  1.3:  Diagram  of  Smoke  Visualization  (adapted  from  [11]) 

locity  at  a  point  requires  calculation  of  the  influence  of  all  vortex  filaments  on  that 
point.  Vortex  methods  have  been  used  in  conjunction  with  other  methods  providing 
inflow  values  to  potential,  Euler  or  Navier-Stokes  solvers. 


1.2.2  Euler  and  Navier-Stokes  Fixed  Grid  Methods 

Conceptually,  Euler  and  Navier-Stokes  solvers  greatly  simplify  the  process  of  calcu¬ 
lating  a  rotor  flow.  They  can  predict  both  the  load  on  the  rotor  and  the  associated 
wake  in  a  uniform  manner  free  of  artificial  constructs  such  as  wake  markers.  In 
actuality,  the  efforts  necessary  to  achieve  adequate  grid  resolution  near  the  wake 
(adaptive  grid  refinement,  overset  and/or  unstructured  grids)  tend  to  complicate  this 
approach.  The  wake  has  been  found  e_xperimentally  to  be  as  thin  as  5%  of  the  chord. 
Supposing  it  requires  12  points  across  the  thinnest  part  of  the  wake  to  adequately 
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1.2.  REVIEW  OF  EXISTING  METHODS 


resolve  the  wake  and  it  is  desired  to  maintain  this  resolution  down  to  one  radius 
below  the  plane  of  a  rotor  with  blade  aspect  ratio  of  15,  then  this  grid  will  require 
2*ir  *  (l5  *  =  2.9  *  1011  points.  And  even  then  the  “solution”  obtained  on  this 

grid  will  be  only  a  model  of  the  Navier-Stokes  equations  for  the  vortex  core  region 
which  can  be  turbulent  with  even  smaller  length  scales.  At  present,  a  rotor  wake 
computation  of  10  million  (1.0  *  107)  points  is  considered  very  large  and  requires  tens 
or  even  hundreds  of  hours  on  expensive,  specialized  supercomputers.  Therefore,  the 
usefulness  of  Euler  and  Navier-Stokes  solvers  is  limited  to  regions  in  the  flow  field  ei¬ 
ther  very  near  the  blade  where  adequate  grid  resolution  can  be  achieved  or  in  regions 
far  enough  removed  from  the  rotor  that  the  rotor  wake  no  longer  has  a  first-order 
impact  on  the  blade  loadings.  The  former  justification  was  used  in  employing  an 
overset  Navier-Stokes  solver  in  the  Helix-i  code  to  determine  the  flow  near  the  rotor 
blade  [12].  The  latter  justification  is  used  in  employing  an  Euler  solver  to  determine 
flow  near  the  ground  in  the  Helix-iii  code. 

1.2.3  Point-Lattice  Methods 

1.2.3.1  Cloud-in-Cell 

The  cloud-in-cell  or  vortex-in-cell  method  utilizes  a  Lagrangian  wake  convected  by 
its  own  induced  velocity  which  is  calculated  on  an  Eulerian  grid.  By  calculating 
the  wake-induced  velocity  on  an  Eulerian  grid  and  interpolating  the  velocities  back 
to  the  wake  in  order  to  convect  it,  an  artificial  viscosity  is  inherently  introduced 
by  truncation  error  which  has  the  benefit  of  eliminating  the  singularity  present  in 
Biot-Savart  methods.  This  is  a  two  dimensional  method.  [13] 
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1.3.  CURRENT  METHOD 


1.2. 3. 2  Vorticity  Embedding 

Vorticity  Embedding  (VE)  is  a  Lagrangian  technique  for  modeling  the  wake  created 
by  a  lifting  body.  Like  other  Lagrangian  methods  it  avoids  the  numerical  diffusion 
associated  with  standard  Eulerian  methods.  VE  satisfies  the  Euler  equations  (in¬ 
tegrated  though  the  wake  sheet)  rather  than  satisfying  the  Biot-Savart  law,  so  its 
usefulness  is  not  limited  to  incompressible  flow  like  traditional  Lagrangian  methods. 
Further,  calculation  of  velocity  at  a  point  does  not  require  calculating  the  influence  of 
every  filament/point/blob  on  that  point.  Therefore,  it  is  computationally  inexpensive 
compared  to  more  traditional  vortex  methods.  [19] 

1.3  Current  Method 

This  particular  method  utilizes  a  lifting  line  together  with  a  VE  solver  to  model  the 
rotor  wake.  This  solver  is  then  coupled  to  a  fractional  step  (FS)  solver  to  model  the 
interaction  of  the  rotor  wake  with  a  ground  plane. 

Helix-III  was  originally  conceived  by  Prof.  Steinhoff  and  implemented  -  without 
the  lifting  line  -  by  Dr.  M.  Moulton  and  Dr.  S.  Babu.  It  was  passed  on  to  the  author 
after  Dr.  Babu  s  departure  from  UTSI.  The  code  required  a  couple  of  corrections  and 
implementation  of  the  lifting  line.  The  required  corrections  and  the  addition  of  the 
lifting  line  are  the  goals  of  this  work. 


Chapter  2 


Numerical  Method 

2.1  Overview 

The  computational  region  is  spanned  by  two  Cartesian  Eulerian  grids  which  have 
one  or  more  coincident  planes  (perpendicular  to  the  z-axis)  as  shown  in  Figure  2.1. 
In  the  “upper”  or  “rotor”  region,  we  utilize  a  lifting  line  to  represent  the  rotor  and 
Vorticity  Embedding  to  model  the  wake  while  in  the  “lower”  or  “ground”  region,  a 
conventional  fractional  step  Euler  solver  is  employed.  Hereafter,  these  regions  will  be 
referred  to  as  the  Vorticity  Embedding  or  VE  region  and  the  fractional  step  or  FS 
region,  respectively.  The  overall  solution  procedure  for  the  Helix-iii  is  outlined  in 
Algorithm  1. 

The  previously  mentioned  Helix-I  code  has  been  used  extensively  for  the  com¬ 
putation  of  hover  flows  where  a  blade-fixed  (rotating)  cylindrical  H-grid  proved  to 
be  most  convenient.  In  a. blade-fixed  grid  (and  neglecting  the  fuselage),  hover  flow 
is  steady  and  blade-to-blade  periodic  making  it  possible  to  compute  the  entirety  of 
the  desired  flow-field  based  on  a  single  blade.  In  forward  flight  or  in  the  presence  of 
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2.2.  LIFTING  LINE  ROTOR  REPRESENTATION 


Algorithm  1  Helix-III  Solution  Procedure  '  " — 

1.  Advance  the  rotor  blade  (lifting  line)  by  A 6  (usually  one  or  two  degrees). 

2.  Perform  the  lifting  line  calculation  to  get  the  circulation  (or  enforce  a  fixed 
circulation)  to  impose  on  the  new  wake  markers  which  are  added  at  the  new 
rotor  position. 

3.  Calculate  the  velocities  in  the  VE  region  using  the  velocity  of  the  FS  region 
from  the  previous  time  step  to  calculate  the  flux  BC  for  the  interface  boundary. 

4.  Calculate  the  velocities  in  the  FS  region  using  the  velocity  from  the  VE  region 
for  the  velocity  BC  at  the  interface. 

5.  Repeat  the  above  steps  until  sufficiently  converged 


a  cross-wind,  such  simplifications  are  no  longer  possible.  Therefore,  Helix-iii  has 
been  devised  based  on  a  uniform  Cartesian  grid.  This  erases  the  distinction  between 
hover  and  forward  flight.  Further,  it  allows  HELIX-III  to  be  coupled  to  another  solver 
in  order  to  better  resolve  flow  near  the  ground. 


2.2  Lifting  Line  Rotor  Representation 


Since  the  grid  is  non-rotating,  the  rotor  must  move  through  the  grid.  Typically, 
this  is  accomplished  with  some  type  of  oversetting  wherein  an  inner  blade-fixed  grid 
communicates  with  the  Cartesian  grid  by  interpolation.  At  present,  this  type  of 
method  is  not  employed  because  only  the  general  flow  field,  especially  near  the  ground, 
is  required.  Therefore,  a  lifting  line  rotor  representation  is  utilized.  In  this  approach, 
the  rotor  blade  is  represented  as  a  line  vortex  and  requires  no  special  gridding.  The 
strength  of  the  sheet  is  the  local  rotor  circulation.  This  is  obtained  by  computing  a 
local  rotor  inflow  angle  (requiring  an  interpolation  from  the  Cartesian  velocity  field 
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2.2.  LIFTING  LINE  ROTOR  REPRESENTATION 


Lifting-line  rotor 
representation 


to  the  lifting  line)  and  assuming  a  local,  ideal  two  dimensional  flow  to  obtain  the 
circulation. 

The  local  angle  of  attack  is  the  sum  of  the  collective  pitch  and  the  inflow  angle 
which  is  calculated  as  follows: 


ai  =  —  arctan 


(2.1) 


where  %  is  the  downward  component  of  the  local  free-stream  velocity  which  may  be 
adjusted  to  include  a  downward  contribution  from  the  bound  vortex  associated  with 
the  leading  edge  of  the  wake  sheet.  That  is, 


=  W‘  +  2rf 


(2.2) 
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2.2.  LIFTING  LINE  ROTOR  REPRESENTATION 


Figure  2.2:  Lifting  Line 


where  d  is  the  distance  to  the  leading  edge  of  the  wake  sheet.  This  requires  a  tri-linear 
interpolation  of  the  downward  component  of  velocity,  wit  at  a  point  one  chord  ahead 
of  the  leading  edge  of  the  wake  sheet.  The  sectional  lift  coefficient  is  then  computed 
according  to  thin  airfoil  theory, 

ct  =  27ra  (2.3) 


where 

=  0.75  +  a*  •  (2.4) 

The  circulation  associated  with  a  marker  originating  at  the  leading  edge  of  the  wake 
sheet  is 

r  =  §  (£»)  •  (2.5) 

This  lifting  line  approximation  provides  reasonable  results  to  within  one  chord  of  the 
blade  tip  where  the  two  dimensional  assumption  is  no  longer  valid.  Therefore,  within 
one  chord  of  the  tip,  the  circulation  is  continued  quadratically  to  zero. 
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2.3.  VORTICITY  EMBEDDING  REGION 


2.3  Vorticity  Embedding  Region 

In  the  Vorticity  Embedding  region,  the  flow  is  determined  by  the  rotor  wake.  It  is 
necessary  to  adequately  resolve  the  rotor  wake  for  at  least  one  blade  passage,  as  there 
is  a  pronounced  effect  on  the  blade  loading.  Therefore,  Vorticity  Embedding  is  used 
in  this  region.  Vorticity  Embedding  is  a  computationally  efficient  way  of  embedding  a 
thin  layer  of  vorticity  (the  wake)  around  Lagrangian  markers  in  an  otherwise  potential 
flow.  This  method  is  the  basis  of  the  Helix-i  code  which  has  been  used  successfully 
for  the  prediction  of  hover  performance  of  isolated  rotors  [15], (16].  The  Lagrangian 
nature  of  the  wake  convection  eliminates  the  need  for  extremely  high  grid  densities 
(to  prevent  diffusion)  and  allows  thin  regions  of  high  vorticity  to  be  convected  over 
distances  large  enough  to  accurately  model  several  revolutions  of  the  rotor  wake. 
The  overall  solution  procedure  for  the  VE  region  is  found  in  Algorithm  2. 

2.3.1  Governing  Equations 

Typical  rotor  tip  speeds  place  the  embedding  region  in  the  compressible  flow  regime. 
However,  for  simplicity  and  to  permit  the  use  of  available  direct  solvers,  the  present 
implementation  is  incompressible.  Therefore,  the  incompressible  version  of  the  con¬ 
tinuity  equation, 

V  •  q  =  0,  (2.7) 

is  satisfied  in  the  VE  region.  The  circulation  introduced  into  the  flow  by  the  lifting 
lines  is  preserved  according  to  Kelvin’s  theorem: 
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Algorithm  2  Vorticity  Embedding _ _ _ _ 

1.  Calculate  the  Clebsch  functions  on  the  grid  nodes  within  the  spreading  distance 
of  the  Lagrangian-marker-wake-sheet. 

(a)  Calculate  the  sheet  strength  (circulation),  V. 

(b)  Calculate  the  shape  function,  A. 

2.  Calculate  the  rotational  component  of  velocity,  qv,  according  to  equation  2.12 

3.  Calculate  the  irrotational  component  of  the  velocity  field,  V <f>. 

(a)  Calculate  4>  such  that 

VV  =  -V  •  qu.  (2.6) 

(b)  The  irrotational  component  is  V<£. 

Therefore,  velocity  (2.11)  satisfies  the  continuity  equation  (2.7). 

(c)  Set  qv  =  qv  +  V<£ 

(d)  Extrapolate  new  boundary  conditions 

(e)  Repeat  (a)-(d)  until  q”  is  sufficiently  divergence-free. 

Then  q  =  qv. 

4.  Stop  if  sufficiently  converged;  otherwise,  convect  the  wake  markers  using  the 
two-step  Runge-Kutta  method  and  start  over. 
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The  circulation  around  a  material  loop  changes  only  at  the  lifting  line,  as  that  loop 
is  convected  through  the  VE  region. 

2.3.2  Vorticity  Embedding 

2.3. 2.1  Complex-Lamellar  Velocity  Decomposition 

It  is  typical  of  velocity  decompositions  to  represent  the  velocity  as  the  sum  of  a 
vortical  component  and  the  gradient  of  a  scalar,  that  is, 

q  =  qu  +  V<£.  (2.9) 

Because  the  curl  of  the  gradient  of  a  scalar  is  zero,  the  vorticity  is  entirely  contained 
within  the  the  vortical  component,  so 

w  =  Vxq=Vxq".  (2.10) 

One  such  decomposition  is  the  complex-lamellar1’2  decomposition.  Here,  the  velocity 
decomposition  takes  the  following  form: 


q  —  rVA  + 


(2.11) 


l<l"  is  a  complex-lamellar  flow  field  if  it  can  be  divided  by  an  integrating  function,  T,  to  produce 
a  potential  or  “lamellar”  flow  [14] : 


2Equivalently,  a  complex-lamellar  vector  field  is  any  vector  field  which  is  perpendicular  to  its  own 
curl  [14],  that  is 

q*-(Vxq*X=0 
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where  the  three  variables  (r,  V,  A)  are  sometimes  called  Clebsch  variables.  The 
complex-lamellar  decomposition  is  very  interesting  in  that  the  Clebsch  variables  can 
be  chosen  such  that  qv  =  0  when  u>  =  0  [14]. 


2.3. 2.2  The  Clebsch  Variables  of  Vorticity  Embedding 

Vorticity  Embedding  is  based  on  a  particular  complex-lamellar  velocity  decomposi¬ 
tion.  The  Clebsch  variables  are  chosen  such  that  the  vortical  component,  q”,  is  only 
non-zero  in  a  thin  region  surrounding  the  Lagrangian  wake  markers.  Said  vortical 
component  is  defined  as 

q«  =  _rVA  (2.12) 

where  T  is  the  circulation  and  A  is  an  as  yet  arbitrary  shape  function.  For  compu¬ 
tational  simplicity,  qv  is  chosen  such  that  it  is  perpendicular  to  the  wake  sheet  [19]. 
Therefore  the  circulation  is 


r  =  Jq’-dn 


(2.13) 


where  the  integration  is  performed  normal  to  the  sheet.  In  terms  of  a  traditional 
potential  flow  solver,  the  integral  of  qv  must  equal  the  local  potential  jump,  T,  across 
the  sheet  of  shed  vorticity.  For  convenience,  A  is  chosen  to  be  a  half-sine  wave  whose 
zero-crossing  is  at  the  sheet: 


(2.14) 


where  Sn  is  the  signed  normal  distance  from  the  sheet  and  a  is  the  spreading  param¬ 
eter  which  determines  the  number  of  cells  over  which  the  vortex  sheet  is  spread  or 
“smeared”. 
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2.3. 2.3  Spreading 

2.3. 2.3.1  Computational  vs.  Physical  Spreading  In  its  original  formulation, 
VE  utilized  physical  spreading.  The  Clebsch  variables  were  calculated  within  a  given 
physical  distance  of  the  wake  sheet.  Using  physical  spreading,  a  search  is  performed 
from  each  Eulerian  grid  point  to  find  the  Lagrangian  wake  markers  that  fall  within 
the  physical  spreading  distance.  The  Clebsch  variables  on  that  Eulerian  grid  point 
are  then  calculated  based  on  a  weighted  average  of  the  contributions  from  all  of  the 
markers  that  fall  within  the  physical  spreading  distance  of  the  grid  point.  The  use  of 
computational  spreading  was  previously  proposed  [17]  and  implemented  [12]  wherein 
the  Clebsch  variables  were  calculated  in  computational  space  rather  than  physical 
space.  Now,  a  represents  a  computational  distance,  or  number  of  cells,  from  the  sheet 
over  which  the  Clebsch  variables  are  calculated.  Additionally,  instead  of  iterating 
over  Eulerian  grid  points,  the  procedure  was  changed  to  iterate  over  the  Lagrangian 
wake  markers  (a  2-D  structure  rather  than  a  3-D  structure). 


2.3. 2.3.2  Weighted  Averaging  Because  each  Eulerian  grid  point  can  be  influ¬ 
enced  by  multiple  wake  markers,  T  and  Sn  (for  the  calculation  of  A  per  2.14)  must  be 
be  treated  as  a  weighted  average  from  all  wake  markers  that  are  within  the  selected 
region  of  influence.  Therefore, 


and 


(Rn)ij,k  — 


E  E44 

.  (  VJt*  ) 


(2.15) 

(2.16) 
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where  i,  j,  k  axe  the  Eulerian  grid  indices,  t  is  the  wake  marker  index,  afj  k  is  the  ratio 
of  the  difference  between  the  spreading  distance  and  the  distance  to  the  grid  point  to 
the  spreading  distance  squared  defined  by 


- max 


(2.17) 


where  is  the  position  of  the  grid  point  and  x*  is  the  position  of  the  wake  marker 
and  A  is  the  sum  of  these  ratios, 

^  =  E  •  (2-i8) 

t  W>*  / 

The  current  numerical  procedure  for  normalization  (multiplying  ritjtk  and  (Sn)ij,k 
by  i)  requires  us  to  iterate  over  all  the  grid  points.  This  defeats  the  intentions 
of  the  change  in  the  Clebsch  procedure  from  iterating  over  grid  points  to  iterating 
over  wake  markers.  The  Clebsch  calculation  again  necessitates  iterating  over  a  3-D 
structure  rather  than  just  over  a  2-D  structure  (the  wake  sheet)  as  might  be  expected. 
Further,  it  has  unfortunate  consequences  when  we  increase  the  number  of  wake  sheets 
to  eliminate  unintended  averaging. 

2.3.2.4  Poisson  Equation 

Given  the  rotational  component  of  velocity,  qv,  we  can  now  calculate  the  irrotational 
component  of  the  velocity,  V<£.  We  require  that  exactly  cancel  out  the  divergence 
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created  by  the  rotational  component  of  velocity.  Therefore, 

V  •  V<£  =  VV  =  -V  •  qu  (2.19) 

providing  us  with  a  Poisson  equation  to  solve  for  the  potential  <f>.  The  regular  grid 
allows  us  to  use  a  fast  direct  Poisson  solver  from  FISHPACK  [see  Appendix  C].  The 
only  other  things  necessary  to  calculate  the  irrotational  velocity  component  are  the 
appropriate  boundary  conditions. 

2.3.2.5  Boundary  Conditions 

In  the  simplest  case:  an  isolated  rotor  in  hover  without  coupling  to  a  ground  region 
(the  FS  region);  Dirichlet  boundary  conditions  are  employed.  The  flow  is  assumed  to 
be  normal  to  the  boundary  surface.  This  is  accomplished  by  using  a  non-homogeneous 
lower  boundary  condition. 

<f>  =  J  qvTdr  (2.20) 

where  q J  is  the  radial  component  of  q°and  the  integration  is  performed  radially,  in 
from  the  outer  edge  of  the  wake  region. 

The  aforementioned  non-homogeneous  lower  boundary  condition  is  no  longer  ap¬ 
propriate  given  the  presence  of  a  ground  (FS)  region  to  which  it  is  coupled.  Instead, 
the  flux  across  the  boundary  is  calculated,  based  on  the  velocity  from  the  coincident 
plane  of  the  FS  region  (from  the  previous  time  step)  and  the  vortical  velocity  from 
the  wake  sheet  according  to  the  following  equation: 

3^  =  (Qfs1  ”  <&£:)  •  n,  (2.21) 
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which  follows  from  the  definition  of  the  total  velocity  in  the  VE  region  as  the  sum 
of  an  irrotational  and  rotational  velocity.  In  this  case,  the  irrotational  component  is 
calculated  from  the  total  velocity  (extrapolated  from  the  interior  of  the  FS  region) 
and  the  rotational  component  from  the  VE  region.  For  all  of  the  other  boundaries,  a 
homogeneous  Dirichlet  conditions  are  initially  imposed  followed  by  repeated  solution 
of  the  Poisson  equation  with  Dirichlet  boundary  conditions  extrapolated  from  the 
previous  solution  until  convergence  is  reached. 


2.3.2.6  Staggered  Computational  Grid 


A  staggered  128x128x32  regular  Cartesian  grid  is  used.  The  scalar  variables  (<f>,  V,  A, 
etc.)  are  stored  at  the  node  points  and  the  vector  quantities  (q  and  qv)  are  located  at 
the  cell  centers.  This  is  convenient  as  gradients  and  divergences  can  be  approximated 
using  a  box  scheme.  In  2-D  for  example. 


Ax 


Ay 


and 


(V*)u  — 


Ax 


Ay 


where  for  example  =  §  (q?+r  M  +  q”+J  j.,)  «"<•  ^  °th«  ^  wilh  tete*“1 
indices  are  similarly  defined. 

In  previous  implementations  of  VE,  the  calculation  of  the  Clebsch  variables  re¬ 
quired  a  computationally  expensive  search  to  locate  the  Lagrangian  wake  points  [12]. 
In  the  present  implementation,  this  search  is  reduced  to  a  simple  mapping  by  the  use 
of  a  regular  Cartesian  grid.  Also,  the  regular  Cartesian  grid  simplifies  the  process  of 
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Figure  2.3:  qv 

calculating  the  Clebsch  variables  as  physical  spreading  and  computational  spreading 
become  synonymous. 

2.3. 2. 7  Two  Dimensional  Implementation 

As  part  of  this  investigation,  a  two  dimensional  version  of  Vorticity  Embedding  was 
briefly  implemented.  Some  results  from  this  implementation  are  presented  here  to 
help  illustrate  the  VE  solution  process. 

Figure  2.3  shows  a  two  dimensional  qu  field  generated  by  a  “lifting  point”  and 
a  trailing  straight  line  wake  which  is  analogous  to  the  branch  cut  in  a  traditional 
potential  flow  calculation.  A  line  integral  of  the  component  of  qu  perpendicular 
to  the  line  of  markers  taken  about  the  edge  of  the  line  of  markers  will  recover  the 
circulation  T  imposed  at  the  leading  edge  of  wake  markers  and  hence  on  all  of  the 
markers  shown.  Figure  2.4  is  the  divergence-free  velocity  field  associated  with  the  qv 
distribution  in  Figure  2.3  and  Figure  2.5  is  that  velocity  plus  a  constant  free-stream 
velocity. 
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Figure  2.4:  qu  +  V<f> 


Figure  2.5:  q 
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2.4  Fractional  Step  Region 

The  Helix-iii  code  will  be  used  for  calculations  in  which  a  ground  plane  is  present.  In 
the  presence  of  a  ground  plane,  the  wake  structure  becomes  complex.  Physically,  we 
see  that  without  a  free  stream  velocity,  the  rotor  wake  system  expands  outward  along 
the  ground  in  rings  of  increasing  diameter.  In  the  presence  of  a  free-stream  velocity, 
the  ground  vortex  forms  a  horseshoe  shape  with  the  ends  of  the  horseshoe  pointing 
downstream.  The  wake  near  the  ground  on  the  free-stream  side  stops  expanding  at 
some  point  as  the  velocity  induced  by  the  wake  becomes  equal  to  the  free  stream 
velocity.  At  this  “stagnation  point”  [6],  the  wake  rolls  up  into  a  concentrated  ground 
vortex.  Lagrangian  methods  (VE  included)  become  cumbersome  when  a  ground  plane 
enters  the  computation.  Not  only  are  wall  boundary  conditions  problematic,  but  a 
computationally  expensive  decision  process  for  combining  the  wake  vortex  elements 
is  necessary.  Fortunately,  the  individual  wake  identities  are  not  important  to  the 
rotor  loads  at  this  distance.  It  was  the  efficiency  of  the  VE  method  for  accurately 
approximating  the  effects  of  these  individual  wake  elements  that  originally  prompted 
the  use  of  VE  in  the  region  near  the  rotor.  Therefore,  a  conventional  solver  is  employed 
in  the  region  near  the  ground. 

2.4.1  Governing  Equations 

In  the  FS  region,  we  assume  that  the  flow  is  incompressible  and  inviscid.  The  equa¬ 
tions  governing  the  physics  of  three-dimensional  incompressible  inviscid  flow  are  the 
incompressible  Euler  equations:  the  momentum  equation, 

9tq  =  ^(q-V)q  — VP,  (2.22) 
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where  P  is  pressure  divided  by  the  constant  density  and  the  continuity  equation, 


V  •  q  =  0. 


(2.23) 


2.4.2  Fractional  Step  Method 

Equations  2.22  and  2.23  (together  with  appropriate  boundary  conditions)  form  a  well- 
posed  set  of  differential  equations  [9];  however,  the  absence  of  an  evolution  equation 
for  the  pressure  makes  constructing  a  straight  forward  numerical  method  difficult. 

The  fractional  step  method 3  (also  commonly  referred  to  as  a  projection  method 4) 
provides  a  way  to  evolve  a  solution  to  the  incompressible  Euler  equations  in  time.  In 
this  case,  the  incompressible  three-dimensional  Euler  equations  are  broken  up  into 
convection  and  mass  balance  equations  which  are  solved  sequentially.  The  overall 
procedure  is  illustrated  in  Algorithm  3. 


2.4.2. 1  Staggered  Computational  Grid 

The  fractional-step  computation  is  also  made  on  a  128x128x32  staggered  grid.  The 
pressure  values  are  stored  at  the  grid  nodes  and  all  of  the  velocity  components  are 
stored  at  the  cell-centers.  A  portion  of  a  two  dimensional  version  of  this  grid  is  shown 

in  Figure  2.6  .  A  similar  grid  was  used  in  [7]  with  a  similar  fractional  step  method  with 

3Temam  uses  this  title.  Temam  proposed  the  method  in  1969.  Yanenko  applies  the  title  to  the 
more  general  class  of  methods  which  includes  ADI,  LOD,  etc.  (21) 

4This  title  was  popularized  by  Chorin’s  justification  for  the  method  as  a  projection  method  in 
the  sense  that  the  velocity  is  projected  from  the  nonsolenoidal  intermediate  velocity  field  into  a 
solenoidal  (divergence-free)  subspace.  [4],  [3]  This  category  of  projection  method  is  not  related 
to  another  category  of  projection  methods  which  includes  Galerkin’s  method  and  the  method  of 
collocation  except  for  their  common  application  to  certain  partial  differential  equations.  Galerkin’s 
method  and  the  method  of  collocation  are  called  projection  methods,  because  they  suppose  a  solution 
exists  in  an  infinite-dimensional  function  space  and  project  the  solution  onto  a  finite-dimensional 
subspace  according  to  basis  functions.  [8j  pg  179] 
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Algorithm  3  Fractional-Step  Method 

1.  The  convection  step 

(a)  Enforce  the  velocity  BC  (taking  the  upper  boundary  velocity  from  the 
rotor  grid) 

(b)  Calculate  the  first  intermediate  velocity,  q* 

q*  =  qn  —  At  (qn  •  V)  qn 

(2.24) 

2.  The  pressure  correction  step  yields  q** 

(a)  Initial  BC  for  P  : 

P  =  0 

and  subsequent  BC  for  P  : 

(2.25) 

P  ~  P zxtrap 

(2.26) 

except  at  the  ground  boundary  where  the  tangency  condition  is  enforced. 
This  leads  to 

DP  1 

9n  =  At  ***  **  ' 

(b)  Calculate  P  such  that 

(2.27) 

V2P  =  -^-V  •  q*. 

At 

(c)  Then  calculate  the  second  intermediate  velocity 

(2.28) 

q"  =  q*  -  At  VP 

(d)  Repeat  (a)  through  (c)  until  P  converges. 

(2.29) 

3.  Apply  smoothing  to  q**,  to  obtain  the  final  velocity: 

qn+1  =  q**  +  cV2q~ 

where  c  is  chosen  arbitrarily  to  keep  the  method  stable  (currently,  c  =  .01). 
The  necessity  of  this  smoothing  step  was  found  by  Dr.  Babu  and  was  accepted 
unquestioningly  by  the  author. 
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Figure  2.6:  The  fractional-step  method  staggered  grid  in  two  dimensions 

the  exception  that  the  pressures  were  located  at  the  cell-centers  and  the  velocities  at 
the  nodes. 

2.4.2.2  Convection  Step 

The  time-discrete  momentum  equation  is 

qn+!  _  qn  _  Af(qn  .  V)qn  -  A tVP"+1  (2-30) 


or 


qn+l  _  q*  _  A tVPn+1 


(2.31) 


where 


q*  =  qn  —  At  (qn  •  V)  qn. 


(2.32) 
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The  calculation  of  q*  constitutes  the  convection  step.  This  calculation  is  performed 
using  a  first  order  “image  point”  method  which  is  attributed  to  (5].  This  method  can 
be  shown  to  be  equivalent  to  the  first  order  explicit  upwind  scheme  for  appropriate 
choice  of  At  (see  Appendix  B  for  explanation  and  equivalence  with  first  order  explicit 
upwind  scheme). 

2A.2.3  Pressure  Correction  Step 

Now,  to  find  Pn+1,  we  take  the  divergence  of  2.30: 

V  •  qn+1  =  V  •  q*  -  AtV2Pn+1  (2.33) 

and  according  to  the  continuity  equation  (2.23),  V  •  qn+1  =  0.  Therefore,  we  are  left 
with  the  pressure  Poisson  equation: 


V2pn+1  =  i_v  .  q.  (2.34) 

Equation  2.34  is  solved  using  a  routine  from  FISHPACK.  The  gradient  of  the  resulting 
pressure  is  then  added  to  the  first  intermediate  velocity,  q*,  in  the  pressure  correction 
step.  The  resulting  velocity,  q**,  then  automatically  satisfies  the  continuity  equation 
(2.23)  as  can  be  seen  here: 


V  •  q**  =  V  •  q*  —  At  (V  •  VP)  =  0.  (2.35) 

Finally,  smoothing  is  added  to  the  second  intermediate  velocity  to  obtain  the  velocity 
at  time  n  +  1. 
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2.4.3  Boundary  Conditions 

The  appropriate  pressure  boundary  conditions  for  the  fractional  step  method  is  a  top¬ 
ic  that  causes  much  confusion  (at  least  for  the  author).  This  is  particularly  true  for 
outflow  boundary  conditions  (see  [18]  for  examples  of  the  variety  of  BCs  employed). 
For  the  outflow  boundaries  (sides,  downstream),  we  approximate  “fully  developed 
flow”,  that  is,  ^  — >  0  by  extrapolating  Dirichlet  boundary  conditions  from  the  inte¬ 
rior.  The  inflow  boundary  conditions  (upstream,  upper)  are  simpler.  The  upstream 
boundary  is  the  specified  free-stream  velocity  (or  is  the  same  as  the  sides  if  there  is 
no  free-stream).  The  upper  boundary  velocities  are  taken  from  the  coincident  plane 
of  the  VE  region.  The  lower  boundary  condition  is  no  flow  through  (tangency)  if  a 
ground  is  present  and  outflow  (as  described  above)  otherwise. 

The  pressure  boundary  conditions  are  also  extrapolated  from  the  interior.  How¬ 
ever,  as  no  evolution  equation  exists  for  the  pressure,  we  are  not  able  to  obtain  an 
interior  solution  without  boundary  conditions.  Therefore,  it  is  necessary  to  assume 
incorrect  homogeneous  Dirichlet  boundary  conditions  in  order  to  obtain  an  interior 
solution.  This  interior  solution  is  then  extrapolated  to  the  boundaries  to  provide 
non-homogeneous  Dirichlet  boundary  conditions.  These  extrapolated  BC  are  then 
presumably  less  incorrect.  This  extrapolation  to  the  boundaries  is  repeated  until 
sufficiently  correct  boundary  conditions  are  obtained.  The  only  exception  to  this  pro¬ 
cedure  is  in  the  case  where  a  solution  with  a  ground  plane  present  is  desired.  In  this 
case,  a  Neumann  boundary  condition  is  found  by  taking  the  dot  product  of  Equation 
2.31  and  the  outward  pointing  normal  yielding 

(2.36) 


d<f> 


^  =  — i-  (qn+1  —  q*V 
dn  At.  \4n  q"/ 


At 
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which  must  be  homogenous  in  order  to  satisfy  the  no-flow-through  condition. 


Chapter  3 


Previous  Results 


One  of  Dr.  Babu’s  original  (unpublished)  results  is  shown  in  Figure  3.1 .  This  result  is 
for  a  Sikorsky  UH-60A  rotor.  This  is  an  untapered  blade  with  an  aspect  ratio,  AR,  of 
15.5.  For  both  cases,  the  tip  Mach  number,  Mt,  is  0.63.  The  blade  circulation  is  fixed 
from  a  previous  solution  using  the  Helix-I  code,  that  is,  the  lifting  line  procedure 
described  in  section  2.2  was  not  used.  The  code  that  generated  the  data  for  these 
images  represents  the  starting  point  for  this  investigation. 


Figure  3.1:  UH-60A  Rotor  in  Hover  with  no  Ground  (Original  Hybrid  Code) 


Chapter  4 


Current  Results 


Two  problems  can  be  seen  in  Figure  3.1.  First,  there  is  a  lack  of  velocity  near  the 
inboard  portion  of  the  wake  which  leads  to  unacceptable  stretching  of  the  wake  sheet. 
And  second,  a  plot  of  the  divergence  reveals  mass  conservation  errors  at  the  interface. 
These  mass  errors  can  also  be  seen  as  unaccounted  for  variations  in  the  velocity  field 
at  the  interface.  As  it  was  also  desired  to  extend  the  code  by  the  addition  of  a  lifting 
line  calculation  to  obtain  the  circulation  to  impose  on  the  leading  edge  of  the  sheet, 
it  was  necessary  to  split  Helix-iii  up  into  two  different  codes:  a  VE-only  code  for 
for  extending  Helix-iii  by  adding  the  lifting  line  and  the  original  hybrid  version  to 
test  changes  to  the  code  to  alleviate  the  two  aforementioned  problems. 

4.1  Vorticity  Embedding- Only  Code 

Vorticity  contours  and  tip  vortex  trajectories  calculated  from  the  results  of  a  hover 
case  run  with  this  VE-only  code  are  shown  in  Figure  4.1.  After  the  results  of  the 
lifting  line  are  examined,  comparisons  will  be  made  between  the  wake  geometry  shown 
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Figure  4.1:  Vortex  Core  Trajectories  and  Vorticity  Contours 
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4.1.  VORTICITY  EMBEDDING-ONLY  CODE 


Figure  4.2:  Variation  of  circulation  due  to  Cartesian  grid 
here  and  wake  geometries  obtained  from  experimental  data. 

4.1.1  Lifting  Line 

4.1. 1.1  Grid  Dependence 

The  interpolation  of  w  introduces  a  small  grid  dependence  in  the  calculated  inflow, 
iWj.  Therefore,  it  was  necessary  to  determine  the  consistency  of  this  approximation 
before  assessing  the  accuracy.  To  this  end,  an  inflow  was  calculated  at  each  of  two 
points  and  recorded  for  one  revolution.  Figure  4.2  illustrates  the  magnitude  of  the 
variation  at  these  two  points  on  the  lifting  line.  The  inflow  calculation  varied  by  less 
than  3%  at  each  of  the  points  and,  therefore,  was  judged  to  be  sufficiently  consistent. 
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4.1.  VORTICITY  EMBEDDING-ONLY  CODE 


Figure  4.3:  Caradonna  &  Tung  Rotor  Circulation,  Helix-III  vs.  Data  (Calculated 
Inflow) 


4.1. 1.2  Comparison  with  Caradonna  &;  Tung  data 

The  circulation,  T,  obtained  by  the  lifting  line  calculation  detailed  in  section  2.2  was 
compared  to  experimental  data  from  [2].  The  results  were  acceptable  for  this  partic¬ 
ular  rotor  (AR=6,  untwisted,  untapered  NACA  0012).  The  lifting  line  assumption 
becomes  invalid  near  both  the  inboard  and  outboard  tips  of  each  blade  and  there 
requires  an  extrapolation  based  on  the  inboard  values  and  the  physical  requirement 
that  the  circulation  goes  to  zero  at  both  the  inboard  and  outboard  tips.  Figure  4.3 
compares  the  circulations  calculated  by  the  lifting  line  in  Helix-iii  with  those  ob¬ 
tained  by  experiment  in  [2].  The  circulation  results  obtained  compare  well  even  for 
this  lower  aspect  ratio  blade  (6). . 
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Circulation 


Figure  4.4:  UH-60A  Circulation,  Helix-I  vs.  Helix-iii  (Fixed  Inflow) 


4.1. 1.3  Comparison  with  HELIX-I  data 

Having  validated  the  lifting  line  representation  of  the  rotor  with  the  experimental  data 
from  [2],  a  lifting  line  computation  was  then  conducted  for  the  UH-60A  rotor.  The 
calculated  circulation  was  then  compared  to  previously  validated  results  from  Helix- 
I.  The  less  than  satisfactory  results  of  this  comparison  are  shown  in  Figure  4:4.  This 
led  to  a.  thorough  re-examihation  and.  re-implementation  of  the  lifting  line  with  the 
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4. 1 .  VORTICITY  EMBEDDING-ONLY  CODE 


same  less  than  satisfactory  results.  This  calculation  was  made  without  allowing  the 
circulation  actually  imposed  at  the  leading  edge  to  vary,  that  is,  the  calculation  was 
made  based  on  an  inflow  that  was  created  by  a  wake  with  the  desired  circulation. 
If  the  circulation  that  is  imposed  at  the  leading  edge  was  changed  according  to  the 
calculated  circulation,  the  circulation  distribution  “evolves”  to  that  shown  in  Figure 
4.5.  The  location  of  the  lifting  line  (one  half  chord  forward  of  the  leading  edge  of  the 
wake  sheet)  was  the  only  parameter  in  the  calculation  that  was  not  taken  strictly  from 
theory,  and  this  was  taken  from  standard  procedures  of  panel  theory  calculations.  It 
was  proposed  that  the  standard  procedures  of  panel  methods  might  not  work  with 
Vorticity  Embedding  and  that  the  appropriate  location  of  the  lifting  line  might  be 
other  than  one  half  chord  in  front  of  the  leading  edge  of  the  sheet.  The  location  was 
iterated  over  from  zero  chords  to  one-and-one-half  chords  forward  from  the  leading 
edge  of  the  wake  sheet  and  up  to  one-and-a-half  chords  above  the  plane  of  rotation. 
The  best  of  these  results  was  still  unsatisfactory.  Attempting  to  account  for  the 
contribution  of  bound  vorticity  to  the  effective  angle  of  attack  according  to  Equation 

2.2  did  not  significantly  improve  the  results.  It  was  then  proposed  that  the  inflow  be 
averaged  azimuthally  about  the  rotor  disk  at  each  radial  station  and  this  inflow  value 
be  used  in  the  lifting  line  computation.  The  slightly  more  encouraging  results  from 
this  further  simplification  are  shown  in  Figure  4.6. 

4.1.2  Comparison  of  Rotor  Plane  Computation  with  Experi¬ 
mental  Results 

The  results  from  the  Helix-iii  code  have  been  compared  to  a  variety  of  data  obtained 
from  experiment.  In  Figure  4.7  ,  the  radial  position  of  the  tip  vortex  computed  by 
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4.1.  VORTICITY  EMBEDDING-ONLY  CODE 


Circulation 


4.1.  VORTICITY  EMBEDDING-ONLY  CODE 


Figure  4.6:  UH-60A  Circulation,  Helix-i  vs.  Helix-III  (Azimuthally  Averaged  In¬ 
flow) 


4.2.  HYBRID  CODE 


Figure  4.7:  Wake  Contraction 

Helix-iii  is  compared  with  the  radial  position  of  the  tip  vortex  measured  in  the 
experiments  of  [10]  and  [2],  A  two-bladed  rotor  with  an  aspect  ratio  of  6  was  used  in 
both  experiments.  The  low  aspect  ratio  makes  the  lifting  line  assumption  invalid  for 
a  large  section  of  the  blade,  but  the  results  still  compare  well  with  experiment.  For 
perspective,  the  contraction  predicted  by  simple  momentum  theory  is  also  plotted. 
In  Figure  4.8  ,  the  vertical  position  of  the  tip  vortex  is  plotted  for  one  revolution. 
The  experimental  data  was  taken  from  the  same  sources  as  above. 

4.2  Hybrid  Code 

In  the  hybrid  code,  the  interface  to  the  FS  region  is  included;  however,  the  lifting 
line  calculation  is  not  attempted.  Instead  a  fixed  circulation,  T  =  r(r),  is  taken  from 
-another  solution.  In  this  case,  the  circulation  .is  taken  from  a  UH-60 A  calculation 
performed  with  the  Helix-i  code.  Two  modifications  were  performed  on  this  code. 
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4.2.  HYBRID  CODE 


Figure  4.8:  Wake  Geometry 


First,  the  wake  was  broken  up  into  segments  to  prevent  overwriting  in  the  calculation 
of  the  Clebsch  variables.  Second,  the  number  of  coincident  planes  between  the  two 
solvers  was  increased  in  order  to  prevent  mass  errors. 

4.2.1  Wake  Segmentation 

The  computational  spreading  approach  to  calculating  the  Clebsch  variables  involves 
iterating  over  each  sheet  of  wake  markers  (a  2D  structure)  rather  than  the  grid  nodes 
(a  3D  structure).  This  results  in  a  significant  savings  in  time,  but  leads  to  complica¬ 
tions  when  resolving  the  wake  for  multiple  revolutions.  Wake  panels  located  further 
down  the  wake  (azimuthally)  wrap  around  and  come  close  to  newer  wake  panels  from 
the  same  sheet.  The  aforementioned  weighted  averaging  (see  Section  2.3.2)  is  no 
longer  appropriate  when  the  wake  panels  are  not  neighboring  panels  but  are  still  part 
of  the  same  sheet.  This  leads  to  the  overwriting  of  T  and  A.  The  effect  of  the  over- 
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Figure  4.9:  Unsegmented  q”  Field 

writing  can  be  seen  as  an  inappropriate  lack  of  velocity  wherever  the  the  sheets  are 
closely  layered  over  top  one  another.  This  necessitates  breaking  up  the  wakes  of  each 
of  the  blades  into  segments  of  less  than  180°.  A  qv  field  generated  by  an  unsegment¬ 
ed  Clebsch  variable  calculation  is  shown  in  Figure  4.9.  A  qv  field  generated  by  the 
Clebsch/qv  routines  after  segmentation  is  shown  in  Figure  4.10.  Breaking  each  sheet 
up  into  segments  and  calculating  the  Clebsch  variables  for  each  segment  rather  than 
each  sheet  ensures  that  non-local  panels  of  a  single  sheet  will  have  an  additive  effect 
on  the  qv  field  rather  than  inappropriately  canceling  each  other  out  according  to  the 
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Figure  4.10:  Segmented  q°  Field 


4.2.  HYBRID  CODE 


averaging  logic. 

An  unexpected  benefit  to  this  wake  segmentation  is  the  ease  with  which  this 
will  allow  for  parallelization  of  the  code.  The  individual  contributions  of  the  wake 
segments  to  the  overall  qv  field  are  independent  of  each  other  allowing  for  simple 
parallelization  for  the  case  where  the  number  of  processors  is  less  than  or  equal  to 
the  number  of  wake  segments. 

4.2.2  Optimization  of  qv  calculation 

The  segmentation  of  each  wake  sheet  leads  to  an  increase  in  computational  time, 
particularly  in  the  calculation  of  qv  from  the  Clebsch  variables.  In  order  to  keep 
computation  times  reasonable  it  was  necessary  to  perform  some  optimization  of  the 
calculation  of  qv.  The  sequence  of  if-then  statements  that  determined  whether  grid 
points  were  above  or  below  the  sheet  or  at  the  edges  was  reduced  to  an  algebraic 
mapping.  The  actual  changes  to  the  code  are  detailed  in  Appendix  D.  These  changes 
reduced  the  overall  computational  time  with  segmentation  by  almost  40%.  Without 
segmentation,  the  reduction  in  computational  time  was  around  15%. 

4.2.3  Boundary  Overlap 

For  the  upper  boundary  of  the  FS  region,  it  is  necessary  to  compensate  for  the  small 
error  introduced  by  the  lower  boundary  condition  in  the  VE  region.  An  overlap  region 
is  created  where  at  least  four  planes  of  the  VE  region  and  the  FS  region  are  coincident. 
The  velocity  boundary  condition  for  the  convection  step  of  the  FS  region  is  then  taken 
from  a  distance  of  five  cells  from  the  lower  boundary  of  the  VE  region.  The  pressure 
boundary  condition  for  the  upper  boundary  of  the  FS  region  is  unaffected- by  the 
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Figure  4.11:  Rotor  in  Ground  Effect.  Calculation  Made  with  Segmentation  but  with- 
out  Increased  Overlap. 

interface  and  proceeds  as  previously  described.  This  change  reduced  the  maximum 
divergence  by  21%.  Qualitative  changes  to  the  velocity  field  are  demonstrated  in 
Figures  4.11  (without  the  increased  overlap)  and  4.12  (with  the  increased  overlap). 

4.2.4  Vortex  Roll-up 

The  eventual  goal  for  this  code  is  to  approximate  the  flow  fields  of  rotors  in  ground' 
effect  in  the  presence  of  a  cross-wind.  To  this  end,  a  calculation  was  made  for  a  UH- 
60A  rotor  in  ground  effect  with  a  cross-wind.  The  results  of  a  particle  trace  through 
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4.2.  HYBRID  CODE 


Figure  4.12:  Rotor  in  Ground  Effect.  Calculation  Made  with  Segmentation  and 
Increased  Overlap. 


4.2.  HYBRID  CODE 


Figure  4.13:  Isolated  UH-60  Rotor  in  Forward  Flight 

this  flow  field  is  shown  in  Figure  4.13.  This  clearly  shows  the  expected  horseshoe¬ 
shaped  ground  vortex. 
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Chapter  5 

Conclusions  and  Recommendations 


Promising  results  have  been  obtained  for  the  hybrid  version  of  the  Helix-III  code. 
The  code  was  compared  to  both  experiment  and  previous  computational  results.  This 
demonstrated  that  the  calculated  wake  geometry  is  acceptable.  The  calculated  lift 
coefficients  were  acceptable  in  one  case  and  less  than  acceptable  in  the  other  case. 
Qualitatively  good  results  were  found  for  the  roll-up  of  the  ground  vortex.  Finally, 
this  paper  will  serve  as  documentation  for  a  code  that  otherwise  may  never  have  been 
documented. 

There  are  a  few  remaining  problems.  First  and  most  important,  it  must  be  possible 
to  calculate  the  correct  circulation  to  impose  on  the  wake  markers.  The  current  lifting 
line  method  is  not  sufficiently  accurate.  The  cause  of  this  is  unknown.  An  overset 
grid  about  the  blade  may  be  a  better  solution.  The  overset  grid  could  take  advantage 
of  the  flow  solver  already  implemented  for  the  FS  region  or  an  established  solver 
could  be  coupled  to  the  embedding  region  as  is  the  case  with  the  Helix-i  code.  To 
simplify  things,  the  could  be  used  only  for  the  calculation  of  circulation,  that  is,  the 
velocities  calculated  in  this  blade  around  the  solver  would  only  be  used  to  calculate 
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the  circulation  to  be  imposed  at  the  leading  edge  of  the  sheet.  Velocity  boundary 
conditions  would  be  interpolated  from  the  Eulerian  grid  and  appropriate  circulations 
would  be  output. 

The  second  remaining  problem  with  the  Helix-III  code  is  speed.  This  lack  of 
speed  is  caused  by  two  things.  First,  the  extrapolation  of  boundary  conditions  for 
the  Poisson  equations  in  both  the  VE  and  FS  regions.  Currently,  five  iterations  of 
the  Poisson  solver  axe  needed  to  extrapolate  values  of  <f>  and  p  at  the  boundaries.  In 
most  fractional  step  methods,  homogeneous  Neumann  BC  are  utilized  for  the  pressure 
calculation.  To  use  Neumann  BC  we  need  the  pressure  defined  at  the  cell-centers  and 
the  velocity  defined  at  the  nodes,  because  FISHPACK  requires  the  forcing  function 
(divergence  of  the  previous  intermediate  velocity)  to  be  supplied  at  the  boundaries 
wherever  Neumann  conditions  are  imposed.  Unfortunately,  because  of  the  coupling  of 
this  solver  to  the  VE  solver,  this  will  require  changes  beyond  the  scope  of  this  work. 

Finally,  the  grid  independence  of  the  solutions  obtained  was  never  demonstrated. 
This  was  pointed  out  towards  the  end  of  the  work  and  was  left  out  due  to  time 
constraints.  However,  this  is  a  basic  step  necessary  for  the  validation  of  CFD  results, 
and  is,  therefore,  the  next  logical  step  for  the  complete  validation  of  the  previously 
obtained  results. 
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Appendix  A 
Nomenclature 


A.l 

Symbols 

AR 

:  blade  aspect  ratio,  * 

a 

:  spreading  distance 

a 

:  angle  of  attack 

Cl 

:  sectional  lift  coefficient 

Mx 

:  hover  tip  Mach  number 

V 

:  pressure 

Q 

:  velocity  vector 

R 

:  rotor  tip  radius 

r 

:  radial  location 

Sn 

:  signed  distance  from  grid  point  to 

wake  marker 

U,  V,  W 

;  velocity  components  in  the  x— ,y— 

and  z —  directions,  respectively 
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Nomenclature 


A.l 

Symbols 

AR 

:  blade  aspect  ratio,  * 

a 

:  spreading  distance 

a 

:  angle  of  attack 

ct 

:  sectional  lift  coefficient 

Mx 

:  hover  tip  Mach  number 

V 

:  pressure 

9 

:  velocity  vector 

R 

:  rotor  tip  radius 

T 

:  radial  location 

Sn 

:  signed  distance  from  grid  point  to 

wake  marker 

u,  u,  w 

;  velocity  components  in  the  x—,y— 

and  z—  directions,  respectively 

52 


Appendix  A 


W{  :  z-component  of  inflow  velocity 
x,  y,  z  :  Cartesian  coordinates 
r  :  bound  circulation 
A  :  shape  function  used  in  VB 
q  :  fluid  velocity  vector,  (u,  v,  w)T 
A  :  uniform  Cartesian  grid  spacing 
A  9  :  rotational  increment 
a  :  spreading  function 
9, 75  :  collective  pitch 
<f>  :  velocity  potential 
X  :  location  of  wake  node 
ft  :  angular  velocity  of  rotating  blade 
to  :  vorticity  vector,  V  x  q 
ip  :  wake  age 

A.2  Superscripts 

*  or  /  :  intermediate  quantity 
oo  :  free-stream  quantity 
l  :  Lagrangian  node 
n  :  time  step 
T  :  rotor  tip  quantity 
v  :  vortical  component 
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Image  Point  Convection 


B.l  Equivalence  with  first  upwind  difference  for  ap¬ 
propriate  At 

The  image  point  method  of  discretizing  the  convection  equation  in  the  first  step  of 
the  fractional  step  method  is  equivalent  to  the  explicit  first  order  upwind  scheme  for 
appropriate  At.  In  one  dimension  on  a  uniformly  spaced  grid,  the  method  proceeds 
as  follows: 

<+i=/(u",xr’*) 


where  the  interpolation  function  I  is  a  function  of  the  velocity  field  at  time  step  n 
and  the  “image”  point  xT**9*  =  X{  —  Aiu"  at  which  u"+1  will  be  calculated.  If  Ui  >  0 
and  At  is  appropriately  small  such  that  |Atu"|  <  Ax  then  I  is  a  function  of  u" 
and  the  image  point. 


Ax 


Ax 


+  E?“- 
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where  rx  =  ^‘.m“se  —  Xi-\  =  xf  —  Atu "  —  a;t-i  =  Ax  —  Atu"1  is  the  distance  between 
the  image  point  and  the  grid  point  to  the  left  of  the  image  point.  Therefore, 


u: 


n+1  _  j  _  A fa"..w 


Ax  W<-1 


Ax  —  Aiu?  _ 


and  after  some  rearranging, 


«r = «?  - 

which  is  the  first  order  explicit  upwind  scheme. 

B.2  Stability 

The  image  point  method  is  a  rare  method  in  that  it  is  both  explicit  and  uncondition¬ 
ally  stable.  It’s  stability  is  derived  from  the  fact  that  the  image  point  at  which  the 
velocity  is  interpolated  always  lies  in  the  appropriate  domain  of  influence  [5].  There¬ 
fore,  it  can  be  used  for  problems  where  the  local  CFL  exceeds  1.0  whereas  the  first 
order  upwind  scheme  becomes  unstable.  However,  if  the  local  CFL  number  exceed- 
s  1.0  near  a  boundary,  the  appropriate  domain  of  influence  (and,  hence,  the  image 
point)  may  lie  outside  the  computational  domain  resulting  in  very  interesting  bus 
errors  or  segmentation  faults.  This  might  be  alleviated  with  “ghost  cells”  outside  of 
the  computational  region,  but  that  point  lies  outside  of  the  theoretical  domain  of  this 
work. 


1Note  that  if  uf  <  0  or  |Atu"|  >_Ax  then  “the  grid  point  to  the  left  of  the  image  point”  will  not 
be  Xi_i  and  the  above  formula  for  rx  will  not  hold. 
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“FISHPACK:  A  package  of 
FORTRAN  subprograms  for  the 
solution  of  separable  elliptic  partial 
differential  equations” 


The  Poisson  equations  (2.6,  2.34)  present  in  both  the  VE  and  FS  grids  are  solved 
using  the  HW3CRT  routine  from  the  FISHPACK  family  of  separable  elliptic  PDE 
solvers  [20].  HW3CRT  solves  the  standard  seven-point  finite  difference  approximation 
to  the  Helmholtz  equation1  on  a  three-dimensional  non-staggered  Cartesian  grid.  It  is 
a  fast  direct  solver  for  the  discrete  equation.  The  solver  first  Fourier  transforms  in  the 
third  dimension  then  utilizes  cyclic  reduction  to  solve  a  reduced  set  of  equations.  A 

non-staggered  grid  is  used,  because  the  scalar  quantities  resulting  from  both  Poisson 

1The  Poisson  equation  is  a  special  case  of  the  Helmholtz  equation:  V2#  +  A<£  =  F{x,  y,  z)  where 
A  =  0.  - 
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equations  are  then  differenced  with  a  box  scheme  which  takes  the  pressure  at  the  grid 
nodes  and  returns  the  gradient  of  the  pressure  at  the  cell  centers  where  the  velocity 
vectors  are  located. 
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qu  Computation  Module 


In  the  qv  routine,  it  is  necessary  to  prevent  large  fluctuations  near  the  edges  of  the 
wake  sheet.  In  order  to  minimize  these  fluctuations,  it  is  necessary  to  determine 
whether  particular  grid  points  lie  above,  below,  or  outside  the  influence  of  the  wake. 
It  is  necessary  to  have  this  information  for  the  calculation  of  the  gradient  of  the  shape 
function,  A,  as  all  of  the  points  involved  in  this  calculation  have  not  necessarily  been 
assigned  a  value  and  will  need  to  be  assigned  a  value  of -.5  or  .5  This  is  accomplished  by 
examining  the  values  of  the  array  facts  (i,  j  ,k)  for  each  of  the  surrounding  points. 
This  array  contains  values  of  -1,  1,  or  0  depending  on  whether  the  grid  point  in 
question  is  located  “above”,  “below”,  or  outside  the  influence  of  the  wake.  To  this  end, 
the  following  lines  of  code  were  executed  at  each  (Eulerian)  grid  point  for  each  wake 
sheet: 

do  141  ii  =1,2 
id  =  ii  +i-l 
do  140  jj  =  1,2 
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jd  =  jj  +j-l 

iw(ii ,  j  j  ,  1)  =  facts  (id,  jd,k) 
iw(ii,jj,2)  =  facts(id,  jd,k+l) 
do  142  kk  =  1,2 
do  139  1  =  1,3 

if (ivv(ii, jj ,kk) .eq.l-2)nl(l)=nl(l)+l 
139  continue 
142  continue 
140  continue 
141  continue 

For  each  grid  point,  nl(l)  is  the  number  of  node  points  of  the  eight  node  points 
surrounding  each  cell  center  which  are  located  “above”  the  wake.  nl(2)  is  the  number 
of  node  points  outside  the  influence  of  the  wake  and  nl(3)  is  the  number  below5  the 
wake.  After  the  introduction  of  segmentation,  this  piece  of  code  was  executed  at  each 
grid  point  for  each  segment  and  required  a  ridiculously  large  number  of  operations. 
This  code  was  replaced  by  the  following  pieces  of  code: 

nl (2)  =8-  (  f acts  (i ,  j , k)  *f acts  (i ,  j  ,k)  + 

&  facts(i+l, j,k)*facts(i+l, j ,k)+ 

&  facts (i+1 , j+1 ,k) *f  acts (i+1 , j+1 ,k) + 

&  facts (i+1 , j+1 ,k+l) *f acts (i+1 , j+1 , k+1) + 

&  facts(i+l, j , k+1) *f acts (i+1, j ,k+l)+ 

&  f acts (i,j , k+1) *f acts (i,j ,k+l)+ 

&  f acts (i , j+1, k+1) *f acts (i, j+1, k+1) + 

ft  f acts (i , j+1 ,k) *f acts (i , j+i ,k)  ) 
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and 


nl  (3)  =  (  (1+f  act s  (i ,  j  ,  k)  )  *f  acts  (i ,  j , k)  + 

&  (1+f acts (i+1 , j , k) ) *f acts (i+1 , j ,k) + 

&  (1+f acts (i+1 , j+1 ,k) ) *f acts (i+1, j+1 ,k)+ 

&  (1+f acts (i+1 , j +1 , k+1) ) *f acts (i+1 , j +1 , k+1)  + 

&  (1+f acts (i+1 , j , k+1) ) *f acts (i+1 , j , k+1) + 

&  (1+f  acts  (i ,  j  ,k+l) )  *f  acts  (i ,  j  ,k+l)  + 

&  (1+f acts (i , j +1 , k+1) ) *f acts (i , j+1 , k+1)  + 

&  (l+facts(i, j+1, k))*f acts (i, j+1, k)  )  /  2 

nl(l)  is  just  8-nl(3)  and  is  no  longer  calculated.  Additionally,  nl(3)  is  only  cal¬ 
culated  if  the  point  in  question  is  not  either  completely  outside  the  influence  of  the 
wake  or  completely  inside  the  influence  of  the  wake  (nl(2)  is  not  equal  to  0  or  8). 
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Jamie  Ryan  Kucab  was  born  eight  days  late  on  April  9th,  1975. 


